library(dplyr)
library(ggplot2)
library(viridis)
library(readr)
library(magrittr)
col_adn <- function() col_factor(c("a", "t", "c", "g", "n"))
col_ADN <- function() col_factor(c("A", "T", "C", "G", "N"))
# 10'' environ
data <- readr::read_csv(
"huge-strong.csv",
col_types = cols(
baseCall = col_adn(),
calledBase = col_ADN(),
uncalledBase = col_ADN(),
## areaUncalledBase = col_double(na = c("-1")),
mutant = col_factor(c("strong", "weak")),
A.base = col_ADN(),
B.base = col_ADN(),
first.call = col_ADN(),
second.call = col_ADN(),
second.area = col_double(),
snp_ = col_factor(c("snp"))),
na = c("-", "", "-1")) %>%
select(
plas = plasmide, # nom de la construction
mutt = mutant, # weak ou strong
base = baseCall, #
ubas = uncalledBase, # base non appelƩe
fcal = first.call,
scal = second.call,
qpos = SCF.position, # position sur la sƩquence
recv = A.base,
donn = B.base,
air1 = area, # aire du pic primaire
air2 = areaUncalledBase,
rel1 = relativeAreaCalledPeak, #
rel2 = relativeAreaUncalledBase,
refp = reference.position,
## upos = positionUncalledBase, # position de la base non appelƩe sur le spectrogramme
## rel2 = parse_double(relativeAreaUncalledBase, na = c("-1", NA)),
## file = file,
## anls = analysis,
## datm = datum,
## fara = first.area,
## sara = second.area,
rat1 = uncorrected.proportion.A,
## rat2 = uncorrected.proportion.B,
snp = snp_,
qual = quality,
comt = comments,
spos = position, # position sur le spectrogramme
## spos1= sposition, # position sur le spectrogramme
A = A,
C = C,
G = G,
T = T
) %>%
mutate(base = toupper(base),
plas = gsub("-1073.ab1", "", plas),
rel2 = parse_double(.$rel2, na = c("-1")),
air2 = parse_double(.$air2, na = c("-1")))
## rƩaffecte les niveaux de facteur.
##
## | facteur | commentaire |
## |---------|-------------------------------------------------------|
## | D | WARNING: bases could not be matched |
## | C | WARNING: multiple peak span, data are not trustworthy |
## | B | WARNING: primary peak does not match the reference |
## | A | processed normally |
## il n'y a pas dans les données de facteur C. d'où les trois niveaux de
## facteur. sensible aux variations de paramĆØtre polySNP Ć mon avis. Ć
## surveiller.
data$comt <- factor(data$comt)
data$comt <- plyr::mapvalues(data$comt,
from = levels(data$comt),
to = c("D", "B", "A"))## # algo :
## par plasmide,
## trouve la première position alignée avec la référence
## ajuste les positions avant celle lĆ en fonction
## ajuste les positions aprĆØs en fonction.
##
## stratƩgie : faire une table qui associe par position la position de rƩfƩrence
## connue la plus proche, inner_join on it, calcul de la diffƩrence entre la
## position de rƩfƩrence la plus proche et la position actuelle
data %>%
filter(!is.na(refp)) %>% # ne garde que les positions ou la ref est connue
mutate(diff = refp - qpos) %>% # calcule la difference entre la position sur
# la sequnece et la position sur la ref
select(plas, diff) %>% # ne garde que les colonnes avec le nom du plasmide et
# la diffƩrence
group_by(plas) %>% # par plasmide
summarise(diff = diff %>% mean() %>% round() ) %>% # compte la moyenne des
# diffĆ©rences, arrondie Ć
# l'entier le plus proche
inner_join(data , ., by = "plas") %>% # regroupe avec la table mĆØre
mutate(apos = qpos + diff) -> # calcule le dƩcalage entre la position sur la
# sƩquence et la position sur la rƩfƩrnece
data
theme_set(
theme_minimal(base_size = 10, base_family = "Courier") +
theme(panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
axis.text = element_text(face = "bold", size = 12)))
## rajoute une colonne longueur qui dƩtermine la longueur de la trace de
## conversion par plasmide.
data %>%
rowwise() %>%
filter(fcal == donn) %>%
ungroup() %>%
group_by(plas, mutt) %>%
summarise(len = max(refp) - min(refp)) %>%
inner_join(x = data, y = .) ->
data
# #' Cette fonction retourne un graphique ggplot2, avec en ordonnƩe les
# #' constructions et en abscisse la position sur la sƩquence de rƩfƩrence.
# #' Chaque point reprƩsente un pic ou phred a dƩtectƩ un pic secondaire. Les
# #' axes verticaux reprƩsentent les positions de snp.
# #' @title plot_second_peak
# #' @param data la table de donnƩe
# #' @param mutt "strong" ou "weak", ne reprƩsente que les positions pour les
# #' constructions weak ou strong
# #' @param length classe les sƩquences en fonction de leur longueur.
# #' @param snp_only ne reprƩsente que les positions de SNP introduits.
# #' @param clean enlève les séquences qui ont trop de polymorphes.
# #' @return un graphique ggplot
plot_second_peak <-
function(data, sw = NULL, length = FALSE, clean = TRUE, snp_only = FALSE)
{
stopifnot(
sw == "weak" || sw == "strong" || is.null(sw),
typeof(length) == "logical",
typeof(snp_only) == "logical"
)
find_length <-
function(data) data %>% arrange(len) %>% {.$plas %>% unique()}
find_snppos <-
function(data) data %>%
filter(!is.na(snp)) %>%
{ if (!is.null(sw)) filter(., mutt == sw) else . } %>%
group_by(refp) %>%
summarise(count = n()) %>% {.$refp}
data %>%
# si clean is true, ne garde que les sƩquences absentes de seqcrade
{ if (clean) filter(., !(plas %in% seqcrade)) else . } %>%
# filter pour ne garder que les sƩquneces strong ou weak
{ if (!is.null(sw)) filter(., mutt == sw) else . } %>%
# filtrer pour ne garder que les snp
{ if (snp_only) filter(., !is.na(snp)) else . } %>%
filter(ubas != "N") %>%
ggplot(
data = .,
aes(x = apos,
# si length, trier par longueur de conversion tract.
y = { if (length) (factor(plas, levels = find_length(data)))
else plas },
color = ubas, fill = ubas, label = ubas)
) +
geom_vline(xintercept = find_snppos(data), color = "gray", alpha = 0.5) +
geom_point(alpha = 0.6) +
xlab("Position sur la sequence de reference") +
ylab("Transformant") +
theme(legend.position = "bottom") +
## scale_color_(discrete = TRUE, name = "Base\nSecondaire") +
scale_color_brewer(name = "Base\nSecondaire", palette = "Set1") +
scale_fill_discrete( name = "Base\nSecondaire")
}
## #' cette fonction dƩfinit si une base est strong ou weak. Si la base en question
## #' est NA, alors renvoit un charactere vide.
## #' @param base : une base d'ADN.
base_caller <- function(base)
{
if (is.na(base) || base == "N" ) { "" }
else if (base == "A" || base == "T") { "W" }
else { "S" }
}
find_length <-
function(data) data %>% arrange(len) %>% {.$plas %>% unique()}
plot_first_second <- function(data, sw = NULL, snp_only = FALSE) {
find_snppos <-
function(data) data %>%
filter(!is.na(snp)) %>%
{ if (!is.null(sw)) filter(., mutt == sw) else . } %>%
group_by(refp) %>%
summarise(count = n()) %>% {.$refp}
data %>%
filter(ubas != "N") %>%
filter(!(plas %in% seqcrade)) %>%
{ if(!is.null(sw)) filter(., mutt == sw) else . } %>%
{ if(snp_only) filter(., snp == "snp") else . } %>%
rowwise() %>%
mutate(type1 = base_caller(base),
type2 = base_caller(ubas)) %>%
ungroup() %>%
ggplot(aes(x = apos,
y = factor(plas, levels = find_length(.)),
label = paste0(type1, type2),
fill = paste0(type1, type2))) +
## geom_tile(alpha = 0.2) +
geom_point(aes(color = factor(paste0(type1, type2))), alpha = 0.5, size = 4) +
geom_text(aes(
color = factor(paste0(type1, type2))),
size = 2,
## alpha = 0.6,
family = "Courier",
fontface = "bold") +
geom_vline(xintercept = find_snppos(data),
color = "gray", alpha = .5
) +
scale_fill_brewer(palette = "Set1") +
scale_color_brewer(palette = "Set1") +
xlab("Position sur la sequence de reference") +
ylab("") +
ggtitle(paste("Constructions", toupper(sw))) +
guides(colour = guide_legend(title = ""), fill = "none") +
theme(legend.position = "bottom",
panel.grid.major.y = element_line(color = "gray", linetype = "dotted"))
}
plot_snp_pos <- function(data, sw) {
faclevel <-
function(x, y) factor(paste0(x, y), levels = c("SS", "SW", "WS", "WW", "S", "W"))
data %>%
filter(snp == "snp") %>%
filter(mutt == sw) %>%
filter(!(plas %in% seqcrade)) %>%
rowwise() %>%
mutate(type1 = base_caller(base),
type2 = base_caller(ubas)) %>%
ungroup() %>%
ggplot(aes(x = apos,
y = factor(plas, levels = find_length(.)),
label = faclevel(type1, type2),
## fill = faclevel(type1, type2),
color = faclevel(type1, type2))) +
geom_point(alpha = 0.4, size = 4) +
geom_text(size = 2, family = "Courier", fontface = "bold") +
## scale_fill_brewer(palette = "Set1") +
scale_color_viridis(discrete = TRUE) +
labs(x = "Position sur la sequence de reference", y = "") +
ggtitle(paste("Transformants", toupper(sw))) +
guides(colour = guide_legend(title = ""), fill = "none") +
theme(legend.position = "bottom",
panel.grid.major.y = element_line(color = "gray", linetype = "dotted"))
}Le graphique suivant reprƩsente toutes les positions par sƩquence auxquelles phred dƩtecte un pic secondaire. La couleur reprƩsente la base au pic secondaire. Les traits verticaux reprƩsentent les positions de SNP introduites. Le trait horizontal dƩmarque les donneurs STRONG (en bas) des WEAK.
plot_second_peak(data, clean = FALSE) +
geom_hline(yintercept = 84.5, linetype = "dotted")Certaines séquences sont extrêmement bruitées. Le graphique suivant représente le nombre de positions polymorphe par plasmide.
data %>%
filter(ubas != "N") %>% # ne regarde que les positions où on a du
# polymorphisme
group_by(plas,mutt) %>%
summarise(ucnt = n()) %T>% # par plasmide, compte le nombre de positions
# polymorphes
{ print( ggplot(data = .) +
geom_point(aes(ucnt, plas, color = mutt)) +
geom_vline(xintercept = 80, color = "red", linetype = "dotted") +
scale_color_brewer(palette = "Set1", guide = FALSE)
)} %>%
filter(ucnt > 80) %>% # au vu de la distribution, 80 semble un bon cutoff de
# propretƩ
{.$plas} ->
seqcrade # sortie dans un vecteur.Jāai utilisĆ© le seuil de 80 positions polymorphes pour discriminer les sĆ©quences propres des sĆ©quences āsalesā (cutoff en rouge). Sans ces mĆŖme sĆ©quences :
plot_second_peak(data, clean = TRUE) +
geom_hline(yintercept = 79.5, linetype = "dotted")## #' Seulement les sƩquences des transformants _strong_
## #+ seqstrong, fig.height = 7.5
## plot_second_peak(data, sw = "strong")
## #' Seulement les sƩquences des transformants _weak_
## #+ seqweak, fig.height = 7.5
## plot_second_peak(data, sw = "weak")Si on ne regarde que les positions de SNP
plot_second_peak(data, snp_only = TRUE) + geom_hline(yintercept = 11.5, linetype = "dotted")On voit plusieurs choses :
data %>%
filter(!is.na(snp)) %>%
filter(!(plas %in% seqcrade)) %>%
filter(ubas != "N") %>%
group_by(mutt) %>%
summarise(count = n()) %>%
ggplot(data = ., aes(x = mutt, y = count, fill = mutt )) +
geom_bar(stat = "identity", width = 0.01) +
geom_point(aes(color = mutt)) +
scale_fill_brewer(palette = "Set1" , guide = FALSE) +
scale_color_brewer(palette = "Set1", guide = FALSE) +
coord_flip() +
labs(x = "Nombre de\npositions polymorphes", y = "Donneurs")data %>%
filter(!is.na(snp)) %>%
filter(!(plas %in% seqcrade)) %>%
filter(ubas != "N") %>%
group_by(mutt, plas) %>%
summarise(count = n()) %>%
ungroup() %>%
arrange(desc(count)) %>%
ggplot(data = .,
aes(y = count,
x = factor(plas, levels = {.$plas}),
fill = mutt, color = mutt)) +
geom_point() +
geom_bar(stat = "identity", width = 0.1, alpha = 0.3) +
coord_flip() +
scale_fill_brewer(palette = "Set1", guide = FALSE) +
scale_color_brewer(palette = "Set1", guide = FALSE) +
labs(x = "Transformants", y = "Nombre de positions polymorphes")data %>%
filter(ubas != "N") %>%
filter(!(plas %in% seqcrade)) %>%
group_by(snp) %>%
summarise(count = n()) %>%
ggplot(data = .,
aes(x = snp, y = count)) +
geom_bar(stat = "identity", width = 0.05, fill = "gray") +
geom_point() +
coord_flip() data %>%
filter(!(plas %in% seqcrade)) %>%
filter(ubas != "N") %>%
group_by(plas) %>%
summarise(totalpoly = n()) %>%
inner_join(x = data, y = . ) %>%
filter(ubas != "N") %>%
filter(!(plas %in% seqcrade)) %>%
group_by(mutt, plas, snp, totalpoly) %>%
summarise(count = n()) %>%
ungroup() %>%
mutate(freq = count / totalpoly) %>%
filter(!is.na(snp)) %>%
arrange(desc(freq)) %>%
ggplot(data = .,
aes(x = factor(plas, levels = {.$plas}),
y = freq)) +
geom_bar(aes(fill = mutt), stat = "identity", width = 0.2) +
geom_point(aes(size = count, color = count)) +
coord_flip() +
scale_color_viridis(begin = 0.2, end = 0.9,
name = "Nombre de positions\navec pic secondaire",
guide = guide_legend(direction = "horizontal",
label.position = "bottom")) +
scale_fill_viridis(discrete = TRUE, begin = 0.8, end = 0.2,
name = "Type de donneur",
guide = guide_legend(direction = "horizontal")) +
scale_size_area(max_size = 10, guide= FALSE) +
labs(x = "Donneur", y = "Frequence des sites avec\npic secondaire correspondant aux marqueurs") +
theme(legend.position = c(0.8, 0.8))data %>%
filter(ubas != "N") %>%
group_by(plas) %>%
summarise(totalpoly = n()) %>%
inner_join(x = data, y = . ) %>%
filter(ubas != "N") %>%
filter(!(plas %in% seqcrade)) %>%
group_by(mutt, plas, snp, totalpoly) %>%
summarise(count = n()) %>%
ungroup() %>%
mutate(freq = count / totalpoly) %>%
arrange(desc(totalpoly)) %>%
mutate(snp = factor(ifelse(is.na(snp), "autre", "snp"),
levels = c("autre", "snp"))) %>%
ggplot(data = ., aes(x = factor(plas, levels= {.$plas}), y = count )) +
geom_bar(aes(fill = snp), stat = "identity", width = 0.3, position = "dodge") +
geom_point(aes(color = snp)) +
scale_color_viridis(begin = 0.2, end = 0.8, discrete = TRUE ) +
scale_fill_viridis(begin = 0.2, end = 0.8, discrete = TRUE ) +
geom_hline(yintercept = 0, color = "gray") +
labs(x = "Donneur", y = "Nombre de positions") +
theme(
axis.text.x = element_text(angle = 90, size = 6),
panel.grid.major.x = element_line(linetype = "dotted", color = "gray")
)Jāai ensuite voulu regarder la rĆ©partition par sĆ©quence de ces seconds pics. La base majoritaire correspond-elle plus souvent Ć la base donneuse, Ć la base receveuse, ou aucun des deux ?
Le graphique suivant reprĆ©sente donc toutes les positions où phred dĆ©tecte un pic secondaire, encore une fois. La couleur du point reprĆ©sente les diffĆ©rents cas de figure : rouge pour les cooccurences SS, bleu pour les cooccurences SW (pic majoritaire S et pic secondaire W), vert pour lāinverse et violet pour les cas WW.
plot_first_second(data)Si on compte le nombre de cooccurences diffƩrentes :
data %>%
filter(ubas != "N") %>%
filter(!(plas %in% seqcrade)) %>%
rowwise() %>%
mutate(altern = paste0(base_caller(base), base_caller(ubas)) %>%
factor(levels = c("SS", "SW", "WS", "WW", "S", "W"))) %>%
group_by(mutt, altern) %>%
summarise(count = n()) %>%
ggplot(data = .,
aes(x = altern, y = count, fill = altern)) +
geom_point(aes(color = altern, size = count)) +
geom_bar(stat = "identity", width = 0.1, alpha = 0.4) +
facet_grid(mutt~.) +
scale_size(guide = FALSE) +
scale_fill_brewer( palette = "Set1", guide = FALSE) +
scale_color_brewer(palette = "Set1", guide = FALSE) +
labs(y = "Nombre de positions polymorphes",
x = "Type de coocurrence") +
theme(legend.position = "bottom") +
coord_flip()Il y a autant de coocurrences WS (donneur / sauvage) quand le donneur est weak que de cooccurrences SW (donneur / sauvage) quand le donneur est strong. Par contre, il y a plus de cooccurrences SW (sauvage / donneur) quand le donneur est weak que de cooccurences WS (sauvage / donneur) quand le donneur est strong. Il y a Ʃgalement un nombre relativement ƩlevƩ de cooccurences SS et WW.
## #+ firstsecondweak, fig.height = 7
## plot_first_second(data, "weak")Si on ne regarde que les positions de SNP pour les donneurs WEAK.
plot_first_second(data, "weak", TRUE)Idem, seulement les SNP pour les donneurs STRONG.
plot_first_second(data, "strong", TRUE)Certains reads ne sont pas homogènes : tantÓt la base donneuse est majoritaire, tantÓt la base receveuse. Le graphique suivant est un alignement de toutes les positions de SNP pour les donneurs WEAK, classés par longueur de trace de conversion.
plot_snp_pos(data, "weak") Idem pour les transformants strong.
plot_snp_pos(data, "strong")On en a conclu plusieurs choses :
il semble que les couleurs associées soit plutÓt bleu avec jaune et violet avec vert. Autrement dit les pics majoritaires dans la trace de conversion sont plutÓt weak chez les donneurs weak, et les pics majoritaires en dehors de la trace correspondent à la base strong.
Dans lāhypothĆØse où ce seraient des contaminations : on dispose du plan de plaque, on peut donc essayer de tester lāassociation entre les contaminations quāon observe et les sĆ©quences prĆ©sentes dans les transformants des autres puits. Une sorte de dĆ©mineur des conta. La mĆ©thode nāest encore pas Ć©tablie.
On retrouve toujours les bases strong (en vert) au milieu des sĆ©quences de weak (en jaune). Ća pourrait ĆŖtre dĆ» Ć des traces de conversions complexes. En tout les cas ces mutations ne sont pas associĆ©es Ć des seconds pics, pour 4 des 6 cas que Laurent avait dĆ©tectĆ© auparavant.
Concernant le problĆØme des contaminations : on sāest dit quāil pouvait subsister lāADN du plasmide, adsorbĆ© sur la bactĆ©rie transformante et quāon isole sur le milieu. La bactĆ©rie isolĆ©e, en cours de croissance, pourrait utiliser cet ADN pour un deuxiĆØme Ć©vĆØnement de recombinaison. Si au lieu dāun Ć©vĆØnement de recombinaison, on a plusieurs Ć©vĆØnements, qui utilisent cet ADN rĆ©siduel comme matrice, Ƨa pourrait expliquer les profils de polymorphismes quāon observe, qui semblent associĆ© (globalement) spĆ©cifiquement Ć nos SNP. On veut donc traiter les clones, avant dāĆ©taler, par de la DNAse, pour Ć©liminer les traces de plasmides rĆ©siduels. # distribution de la longueur de trace de conversion
data %>%
filter(!is.na(snp)) %>% # ne regarde que les positinos de snp
filter(!(plas %in% seqcrade)) %>% # enleve les sequences de faible qualitƩ
rowwise %>% #
filter( base == donn ) %>% # détermine la première base sur la référence qui
# n'est pas une base donneuse.
ungroup() %>%
group_by(plas, mutt) %>%
summarise(switchp = min(refp, na.rm = TRUE) ) %>% # appelƩe switchp
group_by(mutt, switchp) %>%
summarise(count = n()) %>% # dƩtermine la distribution par mutant des ces
# switchp
ggplot(., aes(switchp, count)) +
geom_bar(aes(fill = mutt), stat = "identity",
alpha = 0.5, position = "dodge", width = 0.5) +
geom_point(aes(color = mutt)) +
scale_fill_discrete(guide = FALSE) +
scale_color_brewer(palette = "Set1", name = "Donneur") +
labs(x = "Position de basculement")On ne voit pas une variation affolante entre la longueur des traces de conversions issues des donneurs weak et strong. # Variations des scores de confiance
Chaque fois quāun pic secondaire est dĆ©tectĆ©, un ratio dāaire est calculĆ© entre le pic primaire et le pic secondaire (ratio = 1 => pics Ć©quivalents). Jāai voulu voir la faƧon dont ces ratios Ć©taient distribuĆ©s sur les sĆ©quences. On peut regarder cette variabilitĆ© intra et inter-sĆ©quence.
Le graphique suivant reprƩsente la distribution des scores de confiance par read. La couleur et la taille des points correspondent au ratio des pics (secondaire / primaire, toujours < 1).
data %>%
filter(ubas != "N") %>%
## filter(mutt == "weak") %>%
filter(!(plas %in% seqcrade)) %>%
ggplot(data = .,
aes(x = apos, y = plas)) +
geom_point(aes(size = rel2, color = rel2), alpha = 0.8) +
scale_color_viridis(begin = 1, end = 0, name = "Ratio") +
scale_size_area(max_size = 3, guide = FALSE) +
labs(x = "Position sur la reference", y = "Transformants") +
theme(legend.position = "bottom",
panel.grid.major.y = element_line(color = "gray", linetype = "dotted"))Il semble que les scores varient assez peu au sein dāune mĆŖme sĆ©quence. Cāest ce que le graphique suivant cherche Ć montrer.
data %>%
filter(ubas != "N") %>%
## filter(mutt == "strong") %>%
filter(!(plas %in% seqcrade)) %>%
ggplot(data = ., aes(x = apos, y = rel2)) +
geom_point(aes(size = rel2, color = rel2), alpha = 0.8) +
geom_line(aes(group = plas, color = rel2), alpha = 0.3) +
geom_vline(xintercept = data %>% filter(!is.na(snp)) %>% group_by(refp) %>% summarise(count = n()) %>% {.$refp},
color = "gray", size = 0.1) +
scale_color_viridis(begin = 1, end = 0, guide = FALSE) +
scale_size_area(max_size = 3, guide = FALSE) +
labs(x = "Position sur la reference", y = "Ratio\nPic Secondaire /\nPic Primaire") +
theme(legend.position = "bottom",
panel.grid.major.y = element_line(color = "gray", linetype = "dotted"))Hormis quelques points, les lignes se croisent assez peu, dāune position Ć lāautre. Le graphique suivant reprĆ©sente la mediane et la dĆ©viation absolue Ć la mĆ©diane des ratios de pic, globalement assez ressĆ©rĆ©s autour de 0.15. 0.2 est un score assez conservateur, il correspond clairement Ć un pic secondaire, lorsquāon le vĆ©rifie visuellement sur les spectrogrammes.
data %>%
filter(ubas != "N") %>%
filter(!(plas %in% seqcrade)) %>%
group_by(plas, mutt) %>%
summarise(med = median(rel2, na.rm = TRUE),
mad = mad(rel2, na.rm = TRUE)) %>%
arrange(mad) %>%
ggplot(data = .,
aes(y = plas, x = med, color = med )) +
geom_point() +
geom_errorbarh(aes(xmax = med + mad, xmin = med - mad )) +
scale_color_viridis(begin = 1, end = 0, guide = FALSE) +
theme(legend.position = "bottom") +
labs(x = "MĆ©diane des ratios par read")Jāai voulu regarder sur lāensemble des sĆ©quences, la frĆ©quence des sites avec un pic secondaire par position. LāidĆ©e est de montrer que les pics secondaires sont plus frĆ©quents aux positions de SNP maqueurs quāon a introduit quāaux autres positions.
data %>%
mutate(has_sp = ifelse(ubas != "N", TRUE, FALSE)) %>%
filter(has_sp) %>%
group_by(apos) %>%
summarise(n_ps = n()) %>%
inner_join(x = data, y = .) %>%
group_by(apos, n_ps) %>%
summarise(n_apos = n()) %>%
mutate(freq = n_ps / n_apos) %>%
ggplot(data = ., aes(apos, freq)) +
geom_point(aes(color = freq, size = freq)) +
geom_line(alpha = 0.1) +
scale_x_continuous(breaks =
data %>% filter(!is.na(snp), mutt == "weak") %>%
group_by(refp) %>%
summarise(count = n()) %>% {.$refp} ) +
scale_y_continuous(limits = c(0, 0.30)) +
scale_color_viridis(begin = 0.95, end = 0) +
## scale_y_log10() +
guides(color = "none", alpha = "none", size = "none") +
ggthemes::theme_fivethirtyeight() +
xlab("Position sur la reference") +
ylab("Frequence de pic secondaire") +
ggtitle("Frequence des bases avec pic\nsecondaire a une position donnee") +
theme(
panel.grid.major.y = element_line(linetype = "dotted", color = "gray", size = 0.1),
panel.grid.major.x = element_line(linetype = "dotted", color = "gray", size = 0.4),
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 6),
axis.title.x = element_text(color = "black", size = 8)) Jāai enfin voulu regarder quelles types de transitions Ć©taient favorisĆ©es selon les donneurs.
purpyr_caller <- function(base)
{
if (is.na(base) || base == "N" ) { "N" }
else if (base == "A" || base == "G") { "pur" }
else { "pyr" }
}
## data %>%
## filter(!is.na(snp)) %>%
## filter(base != recv) %>%
## mutate(trans = paste0(base, recv)) %>%
## group_by(mutt, trans) %>%
## summarise(count = n()) %>%
## ggplot(data = .,
## aes(x = trans, y = count )) +
## geom_point() +
## geom_bar(aes(fill = mutt), stat = "identity")
data %>%
filter(!is.na(snp)) %>%
filter(base != recv) %>%
rowwise() %>%
mutate(purpyr_base = purpyr_caller(base),
purpyr_wt = purpyr_caller(recv),
purpyr_trans= paste0(purpyr_wt,"->",purpyr_base)) %>%
filter(purpyr_base != "N") %>%
ungroup() %>%
group_by(mutt, purpyr_trans) %>%
summarise(count = n()) %>%
ggplot(data = .,
aes(x = purpyr_trans, y = count, size = count, color = mutt )) +
geom_point() +
scale_color_viridis(begin = 0.1, end = 0.9, discrete = TRUE) +
scale_size(guide = FALSE) +
labs(x = "", y = "Nombre de transitions") +
coord_flip() +
theme(panel.grid.major.y = element_line(linetype = "dotted", color = "gray"))Le but est de faire un plan de plaque pour Ć©tudier lāassociation entre le nombre de position polymorphe et la position dans la plaque.
plaque <- matrix(1:96, nrow = 8, ncol = 12)
colnames(plaque) <- LETTERS[1:12]
rownames(plaque) <- letters[1:8]
plaque <- reshape2::melt(plaque, value.name = "puit")
colnames(plaque) <- c("colonne", "ligne", "puit")
data %>%
filter(!(plas %in% seqcrade )) %>% # ne garde que les sƩquences propres
filter(!is.na(snp)) %>% # ne garde que les positions de SNP
filter(ubas != "N") %>% # ne garde que les postions de pic secondaire
rowwise() %>%
mutate(plas = gsub(pat= "p.", rep = "", x = plas) ) %>% # remplace les
# caractĆØres inutiles
ungroup() %>%
group_by(plas, mutt) %>% # compte le nombre d'occurence par transformant
summarise(count = n()) %>%
ungroup() %>%
mutate(plas = as.integer(plas)) %>% # mĆŖme types pour le join
left_join(plaque, y = ., by = c("puit" = "plas")) %>% # joint sur le puit
mutate(count = replace(count, is.na(count), 0)) %>% # remplace les NA par 0
filter(!is.na(mutt)) %>% # enlève les positions où on n'a rien
ggplot(aes(y = factor(colonne, levels = c("h", "g", "f", "e", "d", "c", "b", "a")),
x = ligne)) +
geom_point(aes(color = count, size = count)) +
facet_grid(mutt~.) +
scale_color_viridis() +
scale_size_area(max_size = 20) +
guides(
colour = guide_legend(title = "Nombre de positions\navec pic secondaire",
size = 2),
size = "none") +
ggthemes::theme_fivethirtyeight() +
theme(
panel.grid.major.y = element_line(linetype = "dotted", color = "gray"),
panel.grid.major.x = element_line(linetype = "dotted", color = "gray"),
legend.title = element_text(size = 8),
plot.title = element_text(size = 18),
panel.margin = unit(1, "cm"),
strip.text = element_text(size = 14, face = "bold")
) +
ggtitle("Emplacement des transformants\nmontrant des pics secondaires")Il semble quand on regarde vaguement comme Ƨa les sĆ©quences ācontaminĆ©esā, ie avec des pics secondaires, quāelles soient relativement clusterisĆ©es autour de quelques positions. On a donc voulu tester par une simulation si la rĆ©partition est alĆ©atoire ou non. Si elle ne lāest pas, Ƨa peut vouloir indiquer une contamination.
# reconstruit la plaque weak en matrice
data %>%
filter(!(plas %in% seqcrade )) %>% # ne garde que les sƩquences propres
filter(!is.na(snp)) %>% # ne garde que les positions de SNP
filter(mutt == "weak") %>%
filter(ubas != "N") %>% # ne garde que les postions de pic secondaire
rowwise() %>%
mutate(plas = gsub(pat= "p.", rep = "", x = plas) ) %>% # remplace les
# caractĆØres inutiles
ungroup() %>%
group_by(plas) %>% # compte le nombre d'occurence par transformant
summarise(count = n()) %>%
ungroup() %>%
mutate(plas = as.integer(plas)) %>% # mĆŖme types pour le join
left_join(plaque, y = ., by = c("puit" = "plas")) %>% # joint sur le puit
mutate(count = replace(count, is.na(count), 0)) %>% # remplace les NA par 0
{.$count %>% matrix(nrow = 8, ncol = 12)} ->
weak_plaque
print(weak_plaque)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0 1 0 0 0 0 2 0 0 0 0 1
## [2,] 0 0 0 0 0 0 0 0 0 7 0 0
## [3,] 0 0 2 0 17 0 0 0 0 1 15 0
## [4,] 3 0 3 0 4 4 0 1 15 2 1 0
## [5,] 15 1 6 0 1 0 0 0 21 2 9 0
## [6,] 0 0 7 10 8 0 0 0 5 0 0 6
## [7,] 0 1 1 1 3 1 0 0 0 0 0 0
## [8,] 0 0 1 0 0 6 0 0 0 0 0 0
# reconstruit la plaque strong en matrice
data %>%
filter(!(plas %in% seqcrade )) %>% # ne garde que les sƩquences propres
filter(!is.na(snp)) %>% # ne garde que les positions de SNP
filter(mutt == "strong") %>%
filter(ubas != "N") %>% # ne garde que les postions de pic secondaire
rowwise() %>%
mutate(plas = gsub(pat= "p.", rep = "", x = plas) ) %>% # remplace les
# caractĆØres inutiles
ungroup() %>%
group_by(plas) %>% # compte le nombre d'occurence par transformant
summarise(count = n()) %>%
ungroup() %>%
mutate(plas = as.integer(plas)) %>% # mĆŖme types pour le join
left_join(plaque, y = ., by = c("puit" = "plas")) %>% # joint sur le puit
mutate(count = replace(count, is.na(count), 0)) %>% # remplace les NA par 0
{.$count %>% matrix(nrow = 8, ncol = 12)} ->
strong_plaque
print(strong_plaque)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0 0 0 0 0 0 0 0 0 1 0 0
## [2,] 0 0 1 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 2 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 1 0 0 0
## [5,] 0 0 0 0 1 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0
## [7,] 0 0 1 0 0 0 0 0 1 0 0 0
## [8,] 0 1 0 0 1 0 0 0 1 0 0 15
Il semblerait que les sĆ©quences contaminĆ©es soient rĆ©parties uniformĆ©ment dans les sĆ©quences strong, mais quāelles soient clusterisĆ©es autour de deux āsecteurā dans la plaque weak. On a vu avec Laurent pour mettre au point une simulation. Celle ci nous permettrait de voir la rĆ©partition des puits avec voisins contaminĆ©s sous lāhypothĆØse dāune rĆ©partition alĆ©atoire. Si la moyenne du nombre de puits avec voisins contaminĆ©s nāest obtenue quāun petit nombre de fois dans la simulation, alors au risque Ć©gal Ć \(10^{-X}\), avec X le nombre de rĆ©pĆ©titions, on peut dire que la rĆ©partition des puits contaminĆ©s nāest pas dĆ»e au hasard.
library(assertthat) # tests
# une fonction qui génère un sampler entre 1 et n.
random_well <- function(n) sample(1:n, 1)
# une fonction qui initialise la plaque Ć 0 conta par puit
initiate_plate <- function() matrix(rep(0, 96), nrow = 8, ncol = 12)
# une fonction qui transforme l'entrƩe en matrice. par dƩfaut, crƩe une plaque
# vide. doublon. conservƩe par inertie. et flemme.
as.plate <- function(mat = initiate_plate()) matrix(mat, nrow = 8, ncol = 12)
# génère une plaque avec n puits contaminés dedans.
plate_generator <- function(plate = as.plate(), n) {
assert_that(is.matrix(plate), is.numeric(n), n < 96)
# une fonction qui contamine un puit dans la plaque s'il ne l'est pas déjà .
contaminer <- function(mat) {
assert_that(is.matrix(mat))
x <- random_well(8) # tire une ligne au hasard
y <- random_well(12) # tire une colonne au hasard
# si le puit n'est pas déjà contaminé, contamine le.
# sinon, tire un autre puit et contamine le.
if (mat[x, y] != 1) { mat[x, y] <- 1; return(as.matrix(mat)) }
else { contaminer(mat = as.matrix(mat)) }
}
# si la plaque n'a pas déjà n puits contaminés, contamine une plaque avec n-1
# puits, et ainsi de suite.
# si la plaque a déjà n puits contaminés, renvoit là .
if (sum(plate < n )) { plate_generator(as.plate(plate) %>% contaminer(), n - 1) }
else { return(as.plate(plate)) }
}
# gƩnƩre n plaques contaminƩes avec n_conta inside.
n_plate_generator <- function(n_plate, n_conta)
lapply(1:n_plate, function(i) plate_generator(n=n_conta))
### example use :
### =============
## n_plate_generator(3, 3)
# une fonction qui teste si les puits voisins sont contaminƩs. Elle est
# profondƩment testƩe. DO NOT EVEN THINK OF TOUCHING IT.
has_nearby_conta <- function(plate, ligne, colonne) {
assert_that(is.matrix(plate), is.numeric(ligne), is.numeric(colonne))
x <- ligne
y <- colonne
p <- plate
nearby <- c()
bord_gauche <- function(col) ifelse(col == 1 , TRUE, FALSE)
bord_droite <- function(col) ifelse(col == 12, TRUE, FALSE)
bord_sup <- function(row) ifelse(row == 1 , TRUE, FALSE)
bord_inf <- function(row) ifelse(row == 8 , TRUE, FALSE)
coin_sup_gauche <- function(row, col) ifelse(bord_gauche(col) & bord_sup(row), TRUE, FALSE)
coin_inf_gauche <- function(row, col) ifelse(bord_gauche(col) & bord_inf(row), TRUE, FALSE)
coin_inf_droite <- function(row, col) ifelse(bord_droite(col) & bord_inf(row), TRUE, FALSE)
coin_sup_droite <- function(row, col) ifelse(bord_droite(col) & bord_sup(row), TRUE, FALSE)
if (coin_sup_gauche(x, y)) { nearby <- c(nearby, p[x, y+1], p[x+1, y+1], p[x+1, y])}
else if (coin_inf_gauche(x, y)) { nearby <- c(nearby, p[x-1, y], p[x-1, y+1], p[x, y+1])}
else if (coin_sup_droite(x, y)) { nearby <- c(nearby, p[x, y-1], p[x+1, y-1], p[x+1, y])}
else if (coin_inf_droite(x, y)) { nearby <- c(nearby, p[x, y-1], p[x-1, y-1], p[x-1, y])}
else if (bord_gauche(y)) { nearby <- c(nearby, p[x-1, y], p[x-1, y+1], p[x, y+1], p[x+1, y+1], p[x+1, y])}
else if (bord_droite(y)) { nearby <- c(nearby, p[x-1, y], p[x-1, y-1], p[x, y-1], p[x+1, y-1], p[x+1, y])}
else if (bord_sup(x)) { nearby <- c(nearby, p[x, y-1], p[x+1, y-1], p[x+1, y], p[x+1, y+1], p[x, y+1])}
else if (bord_inf(x)) { nearby <- c(nearby, p[x, y-1], p[x-1, y-1], p[x-1, y], p[x-1, y+1], p[x, y+1])}
else {
nearby <- c(nearby,
p[x-1, y-1], p[x, y-1], p[x+1, y-1],
p[x-1, y ], p[x+1, y ],
p[x-1, y+1], p[x, y+1], p[x+1, y+1] )
}
## normalise to 1 if > 1
nearby[nearby > 0] <- 1
## return TRUE if one neighbor is contaminated, false either.
## ifelse(sum(nearby) > 0, TRUE, FALSE)
## compte la moyenne du nombre de cas positif pour ce puit
return(mean(nearby))
}
# compte la moyenne par plaque de la moyenne du nombre de puits DONT LES
# VOISINS sont contaminƩs
count_nearby_conta_well <- function(plate) {
# retrouve les indices des positions par colonne qui sont contaminƩs
locate_conta_well <- function(plate) which(plate != 0, arr.ind = TRUE)
# listes des coordonnƩes des puits avec conta
well_list <- locate_conta_well(plate)
# applique Ć cette liste la fonction qui permet de retrouver la moyenne du
# nombre de puit contaminƩ autour de la position donnƩe par x[1] et x[2]
apply(well_list, 1,
function(x) has_nearby_conta(plate, ligne = x[1], colonne = x[2])) %>%
# moyenne ces moyennes sur l'ensemble de la plaque
mean()
}
# compte le nombre de puits contaminƩs
count_conta <- function(plate) length(which(plate != 0))
count_conta(strong_plaque) # => 11## [1] 11
count_conta(weak_plaque) # => 35## [1] 35
# nombre de plaque dans la simulation
n_plaque <- 10000lance la simulation pour les plaques.
## crƩe une liste de n_plaque qui contient count_conta(*_plaque) rƩparties
## alƩatoirement sur la plaque.
random_weak <- n_plate_generator(n_plaque, count_conta(weak_plaque) ) # takes approx 20''
random_strong <- n_plate_generator(n_plaque, count_conta(strong_plaque)) # takes approx 20''
random_weak_mean <- lapply(random_weak , count_nearby_conta_well) # takes approx 30''
random_strong_mean <- lapply(random_strong, count_nearby_conta_well) # takes approx 30''
# une fonction pour convertir les rƩsultats de lapply en vecteur puis data.frame
# et reprƩsenter la distribution.
plot_mean_distrib <- function(rdm_list) {
assert_that(is.list(rdm_list))
rdm_list %>% unlist() %>% data.frame(wellcount = .) %>%
ggplot(aes(wellcount)) + geom_histogram(binwidth = 0.0005, fill = "gray") +
labs(x = "Moyenne par plaque de la moyenne du\nnombre de puits voisin contamine par puit", y = "") +
theme(panel.ontop = TRUE,
panel.grid.major.y = element_line(colour = "white", size = 0.5, linetype = "dotted"),
panel.grid.minor.y = element_line(colour = "white", size = 0.5, linetype = "dotted"))
}plot_mean_distrib(random_weak_mean) +
geom_vline(xintercept = count_nearby_conta_well(weak_plaque), color = "red") plot_mean_distrib(random_strong_mean) +
geom_vline(xintercept = count_nearby_conta_well(strong_plaque), color = "red") # en admettant une distribution normale, calcule la statistique associƩe
pnorm(q = count_nearby_conta_well(weak_plaque),
mean = mean(random_weak_mean %>% unlist()),
sd = sd(random_weak_mean %>% unlist()),
lower.tail = FALSE) %>% print() # => 0.003313558## [1] 0.003663744
pnorm(q = count_nearby_conta_well(strong_plaque),
mean = mean(random_strong_mean %>% unlist()),
sd = sd(random_strong_mean %>% unlist()),
lower.tail = FALSE) %>% print() # => 0.8288228## [1] 0.8339551